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The equation of state of a system at equilibrium may be derived from the canonical or the grand 
canonical partition function. The former is a function of temperature T, while the latter also 
depends on the chemical potential fi for diffusive equilibrium. In the literature, often the variables 
/3 = (fc^r) 1 and fugacity z — exp(/3/i) are used instead. For real /3 and z, the partition functions 
are always positive, being sums of positive terms. Following Lee, Yang and Fisher, we point out 
that valuable information about the system may be gleaned by examining the zeros of the grand 
partition function in the complex z plane (real /3), or of the canonical partition function in the 
complex f3 plane. In case there is a phase transition, these zeros close in on the real axis in the 
thermodynamic limit. Examples are given from the van der Waal gas, and from the ideal Bose gas, 
where we show that even for a finite system with a small number of particles, the method is useful. 




I. INTRODUCTION 

In a senior level undergraduate or a beginning gradu- 
ate course, examples of phase transition arc often given 
from the classical van der Waals equation of state, and 
at a quantum level, from the Bosc-Einstcin condensation 
(BEC) of an ideal Bose gasi^ The treatment, naturally, 
focuses on the equation of state of the system as a func- 
tion of physical parameters like the temperature T and 
chemical potential //, which are real. It is instructive to 
learn, however, that the approach towards phase transi- 
tion may be studied by examining the analytical behav- 
ior of the partition function for complex values of the 
parameters, even for a finite system where there is no 
discontinuity in the derivatives of the free energy. 

For a system in thermal and diffusive equilibrium, it 
is convenient to calculate the ensemble average using 
the grand canonical partition function Z(Ji, z), where 
ft = {ksT)^ 1 , ks being the Boltzmann constant and 
z = exp(/3^) the fugacity. Note that there is an implicit 
volume dependence in Z, since the eigenenergies are vol- 
ume dependent. We suppress this in our notation Z{(3, z) 
for simplicity. The grand canonical partition function is 
a sum of positive definite terms for real positive values of 
/3 and z, and as such cannot have a zero in the physical 
domain of these variables. Lee and Yang^^ considered a 
lattice gas with a hard core interaction. Because of the 
short-range repulsion between the particles, only a finite 
number of particles may be packed into a finite volume. 
As we shall see, this allows one to express Z as a finite- 
degree polynomial in fugacity z. This polynomial is then 
completely defined in terms of its zeros on the complex 
fugacity plane. These zeros are all complex, coming in 
complex conjugate pairs. In the thermodynamic limit, 



Part of this paper is based on the results reported in the unpub- 
lished undergraduate theses of Calvin Lobo and Allison MacDon- 
ald. 



the zeros coalesce in continuous lines, tending to pinch 
the positive real z axis at a phase transition. In this 
paper, we show that this tendency sets in even at finite 
particle number and volume, with the zeros moving closer 
to the real axis as the particle number is increased. Even 
though complex, the closer a zero comes to the real axis, 
the more it dominates real thermodynamic properties. 
Note that the validity of the Lee- Yang method rests on 
the repulsive core between the particles. 

Fisher— pointed out that the zeros of the canonical par- 
tition function Z?j(f3) on the complex /3 plane have an 
analogous behavior. However, Fisher zeros can give use- 
ful information even in the absence of the repulsive core. 
The Fisher zeros for ideal trapped bosons were studied 
in the context of BEC by Mulken et a/.— In this paper, 
following the work of Hemmer et a/.,— we study the Lee- 
Yang zeros of the classical van der Waals gas (that has 
a phase transition with a critical temperature) and com- 
pare it with a Calogero gas^ (that has no phase tran- 
sition). The ideal Bose gas (which undergoes a phase 
transition at BEC) is studied in some detail, both us- 
ing the grand canonical and the canonical formalisms. 
The heat capacity per particle in the grand canonical 
and canonical ensembles are compared at BEC to check 
how close these are for finite particle number. Since there 
is no short-range repulsion in the ideal Bose gas, the Lee- 
Yang zeros are not meaningful, but the Fisher zeros are. 
Accordingly, we find the pattern of these on the com- 
plex /3 plane for 50 and 100 atoms. Even for such small 
number of particles, we sec a clear tendency for the ze- 
ros to close in on the real /3 axis. The calculations are 
done for the exact Z^{P) as well as the more commonly 
used continuous density of states. This will be discussed 
after the patterns of zeros are presented. The grand par- 
tition function, on the other hand, is shown to have a 
pole at phase transition for real z = L 10 ' 11 A nontriv- 
ial modification over the ideal gas will be to introduce 
interparticlc interaction through virial coefficients in the 
grand potential^ and study how the zeros of the grand 
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canonical partition function shift on the complex plane. 
This is beyond the scope of the present paper. 



II. PARTITION FUNCTION 

The canonical partition function (for fixed JV) is de- 
fined as 



Z N (p) = ]T exp(-/3 J B| Ar) ) , 



(1) 



where E^ N ' are the complete set of eigenenergies of the 
JV-body system including states in the continuum, if any. 
The sum is taken over all states i including the degen- 
eracies. Since the energies depend on the volume of 
the system, there is a volume dependence in Z^(/3). The 
grand canonical partition function, on the other hand, 
allows for particle exchange (in addition to energy) via 
the reservoir, and is defined as 



Z(/3,z)= £ $>xp( 
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N=0 



(2) 

where the fugacity z = exp(/3/z). The simplest example 
is an ideal classical gas in a volume V. The JV-particlc 
canonical partition function is the JVth power of the one- 
particle partition function Z\(f}). The latter is calculated 
by integrating exp(— p 2 /2m) over the phase space divided 
by h 3 , where h is the Planck's constant. The net result 
is 
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where At 
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is the thermal wave length, and 
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the customary division by JV! has been made to preserve 
the extensive property of the entropy. By substituting 
Eq. © in Eq. ©, we obtain 



±UZM(p,z) 



(4) 



where the superscript on Z denotes a non-interacting sys- 
tem. 

More generally, for an interacting gas with short-range 
intcrparticle repulsion, Eq. ^ shows that Z is a finite 
degree polynomial in z, and therefore may be completely 
defined in terms of its zeros. These zeros, however, can- 
not be on the real positive z-axis, since every term in 
Eq. @ is then positive. Accordingly, for complex z (but 
real positive P), Eq. (j4]) may be generalized to^ 



z(p,z)=n(i 



(5) 



The zeros z r (/3, V) come in complex conjugate pairs since 
the coefficients Z N (P) of the polynomial ([2]) are real. In 
case there is a phase transition at some temperature, a 
zero and its complex conjugate tend to pinch the real z 
axis. If there are more than one phase transitions, there 
are segments on the real positive z-axis that are zero- 
free. The grand potential fi = -k B T\nZ(P,z) = -PV, 

hence PP = — \nZ(z,P). Yang and Leo^ proved in gen- 
eral that for V — > oo, this limit exists, and P increases 
monotonically in the zero-free segments. At the inter- 
face of two phases, the pressure P remains continuous, 
but its slope as a function of In z is not the same. More- 
over, in the thermodynamic limit, the number density 
p = g-^j— {PP) may be discontinuous as a function of 
In z in the interface of two phases (see FinkelsteinJ^) 

For our example in the next section, it is more relevant 
to take a finite volume V which contains JV particles. If 
there is a hard-core repulsion, then there is a maximum 
number JV max that can be accommodated in this volume. 
Then the infinite upper limit in the sum over JV in Eq. ([2|) 
is replaced by JV max , and Z(P, z) is a polynomial of order 

JV max • 



A. van der Waals gas 

This is a classical example, studied in the context of 
Yang-Lee zeros by Hcmmcr et aZ.~ Equation ([2|) is used 
with the canonical partition function that is postulated 
to be 



Z N (P) - ±(V 



Nb) exp 



2aP N(N - 1) 
~V 2 



(6) 



where b is interpreted as the volume associated with the 
repulsive core, and a as a measure of the outer attrac- 
tion. For calculating Z{P,z) using Eq. ([2]), one takes 
JVmax = V/b. Note that the excluded volume effect, and 
the outer pair-wise attraction are both incorporated in 
the canonical Z^(P) above. The equation of state can 
easily be deduced from the postulated Zm{P)- One ob- 
tains the Helmholtz free energy F = —\/P\ivZm{P), and 

(dF\ 
— — ) . A little algebra then yields 
dVJ T 

the equation of state 



JV 

P + y^a) {V - Nb) = Nk B T . (7) 



where we have assumed JV 3> 1. One makes the Maxwell 
construction across the unphysical region in which P de- 
creases with V to obtain the equation of state. The criti- 
cal point is obtained by additionally imposing the condi- 
tion that the first and second partial derivatives of P with 
respect to V (at constant temperature ) are zero; see for 
example Landau and Lifshitzii The critical temperature 
is given by 
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FIG. 1. The Yang-Lee zeros of the grand partition function for a van der Waals system of 40 particles, i.e., V = 40, d = 1, 
N — 40. From left to right the patterns correspond to v — 2, v — 27/8, and v = 4. 
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Our interest here is to compute the zeros of Z{(3, z) using 
Eqs. © and © in the complex z-plane for real values 
of T. From Eq. ([6]), we see that there is a cut-off in 
the upper limit of the summation over N = N max = 
V/b, which prevents us from obtaining Z analytically. 
Note, from Eq. ([5]), that a/b has the dimension of energy. 
Setting a = vksl ', and 6=1, Eq. (JSJ takes the form 
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For our calculations, we take -ZV max = V/6 = 40. In 
Fig. [U we plot the zeros of Z(T, z) for three choices of v, 
corresponding to T > T c , T = T c , and T < T c . Even for 
N max = 40, we clearly see the zeros closing in on the real 
z-axis for T < T c . system of 40 particles, i.e., V = 40, 
d = 1, N = 40. From left to right the 



B. Calogero gas 

The Calogero gas is an exactly solvable one- 
dimensional model where point particles are interacting 
with a pair-wise inverse-square potential£i^ The particles 
are trapped in a harmonic oscillator (HO) potential. For 
a repulsive interaction, the high temperature limit of the 
canonical partition function is given byi^ (H = 1) 



1 1 

AH {(3w) N 



exp (- oc/3ujN(N -l)/2) , (10) 



where u> is the oscillator frequency, and a is a measure 
of the strength of the inverse-square two-body poten- 
tial. Since the density of states is a constant in a one- 
dimensional HO, it is like a two-dimensional gas. The 
oscillator length is I = y h/Mui, and we may define a 
density n = N/l 2 = Nu, with h = M = 1. The thermo- 
dynamic limit is taken as N — > oo, uj — > 0, with Nu = n 
a constant. For N ^> 1, we then get 



Zn(P) 
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cxp(-a/?n7V/2) 



(11) 



With this Zn((3), the grand partition function Z(f3,z) 
may be obtained analytically by summing over all N —> 
oo . This is so because these are point particles with no 
excluded volume. A little algebra immediately gives 



lnZ(/3,z) 



/3w 



exp(— na/3/2). 



(12) 



This is of the same form as Eq. © of the perfect classical 
gas, monotonically increasing with z. There is, of course, 
no phase transition. 



III. IDEAL TRAPPED BOSONS AND BEC 

Quantum effects are manifest in a gas when the de 
Broglie thermal wavelength of a particle is larger than, 
or of the order of, the average interparticle spacing. 
BEC was first experimentally realized when neutral 87 Rb 
atoms^ and 23 Na atoms^ were magnetically trapped in 
a HO potential at a few hundred degrees nano-Kelvin. 
For a large number of identical bosons, a sizable frac- 
tion of them abruptly start occupying the lowest level 
even at a temperature much, much larger than the en- 
ergy spacing ttw. This is the condensation temperature. 
We follow the treatment of Ketterle and van Drutenii 
to find the behaviour of the chemical potential /j as the 
temperature is lowered. Here both the temperature and 
chemical potential are taken to be real. Consider N ideal 
bosons in the grand canonical ensemble occupying a dis- 
crete spectrum of single-particle states with energies e n 
at temperature T. Its grand partition function may be 
written as^ 



Z(0,z)=]l{l-zexp(-Pe n ))- 



(13) 



We then get 
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where Z\{0) is the exact one-particle partition function 

and z = e^. Since the occupancy factor of a state has 
to be positive, it follows from the first term on the RHS 
that the smallest value can take is eq, the lowest energy 
single-particle state. In the following, we choose £q = 0, 
so that z < 1, and the power series in z does not involve 
large numbers. Up till now the formulae are general. 
We now specialize to a shifted harmonic oscillator energy 
spectrum, which in one dimension is given by nfius, with n 
going from zero to oo. For an isotropic three-dimensional 

harmonic oscillator z[ 3d} (/3) = (z[ ld) (pj) . Thus we 
write for the three-dimensional HO 



oo / 
- X 1 



OO 

1=1 



(15) 



Note that with the choice of zero-energy ground state, 
z[ ld \/3) = 1/(1 - exp(- / 9w)). In Fig. [2j we plot the 
variation of the chemical potential and the fugacity as a 




Noting that ln(l — x) = — — , a few steps give 



i=i 



]hz(p,z) = J2t Zi ^ 



(18) 



i=i 



where 



Z 1 ^) = (l-expH^))- 



(19) 

The derivative of (E) with respect to T gives the heat 
capacity. 

In Fig. [3j the heat capacity per particle at constant u> 
is shown for (N) = 500 and (N) — 50 based on the grand 
partition analysis of this section. Next we shall compare 
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FIG. 3. The heat capacity per particle for a system of trapped 
bosons as a function of temperature. The calculations based 
on the grand partition function, i.e., Eq. (|16[) , are labelled 
"grand". The graphs labelled "canon" are based on the 
canonical partition function discussed later. 




FIG. 2. The chemical potential fi and the fugacity 2 as a 
function of temperature T when eo ~ for a system of bosons 
in a three-dimensional HO with (N) = 500. 

function of T for a system of trapped bosons in a three 
dimensional isotropic harmonic oscillator when (N) = 
500 (we have set u = 1). Note that the constraint on 
(N) makes fj,, and therefore z temperature dependent. 
There is no discontinuity in /1 or z for a finite number of 
particles, but there is a hint of rapid turning to a plateau 
in both cases near T = 7. This gets more pronounced as 
(N) gets larger. Similarly we can calculate the average 
energy in the grand canonical ensemble, 



<£> = E 



00 ~ l d 



,/3(e« - M) - 



1=1 



For later use, we write the following expression for 
hxZ{p,z) from Eq. (|T3J), 

lnZ(/3,z) = -^ln(l-zexp(-/?e n )) (17) 



these results with the canonical formalism, which requires 
knowledge of Zn(/3). For finite N and (N), the canon- 
ical and grand canonical ensembles may yield different 
results. 

Before concluding this section, we note that for ideal 
bosons, in the thermodynamic limit, there is an analyt- 
ical simple pole at z = 1, the condensation point. This 
may be seen from Eq. (fl~3j). which shows that there are 
poles at 

z„ = exp(/3e„). (20) 

With our choice of e„ > 0, the RHS above is > 1. But 
as the inset of Fig. [2] shows, physically allowed z < 1. 
Therefore from Eq. ([2U|) . the only pole in the physical 
region is at z = 1. 



A. Calculation of Zn{P) 

In order to obtain the grand partition function Z(f3, z) 
from Eq. ((5J), we need to calculate the canonical parti- 
tion function Zm{P)- As we shall see soon, the zeros of 
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FIG. 4. The Fisher zeros for a systems of 50 or 100 trapped bosons taking into account the exact discrete energy spectrum. 



the canonical partition function on the complex (3 plane 
(called the Fisher zeros) are interesting in their own right. 
For an ideal boson or fermion gas, Z^ifi) may be ob- 
tained from Zi((3) using a recursion relation^ We give 
an outline of the derivation of this important relation for 
ideal bosons. We start with the relation (fl8|) . On fur- 
ther expanding the exponential in a power series, and 
equating power by power to the series given by Eq. (|2|), 



Z(/3, z) = 1 + zZ x (fi) + z 2 Z 2 {P) + z 3 Z 3 (P) 



(21) 



the desired recursion relation emerges, which for bosons 
is given by-i^ 



ZnW) 



1 N 

-Y 

n=l 



Zi(n0)Zjv_„(/3). 







(22) 



It is now straightforward to obtain (E) = hi Z^{f3), 

and its derivative with respect to T to obtain the heat 

a 1 d{E) 



capacity. In Fig. [3] we plot = 
labelled "canon". 



as the graphs 



Fisher zeros for ideal bosons 



(the contribution of the state at e = ), F does not 
change sign, and always remains negative. Nevertheless, 
Zjv(/3) is an entire function of (3 in the complex plane. 
For trapped bosons in a 3-dimensional HO, we may use 
the exact Z\(fi) given by Eq. (|19p . Then, using the re- 
cursion relation (|22|) . Zn((3) is a polynomial in the vari- 
able y = in this example of a 3-dimensional HO.~ 
(See Appcndix[A]) The exact form of ^i(/3) is taken, but 
whether Eq = 3/2 or e — the zeros occur at the same 
positions in the complex (3 plane. Only a small fraction 
of the zeros are shown in the Fig. [4] 

Even for N — 50, there is a tendency for the zeros to 
approach the real /3-axis, but clearly the system is not 
condensed. We also display the plot of zeros for N = 100 
ideal bosons. Calculations are much easier if one uses 
Zi(f3) = 1//3 3 , which is the leading term of the exact 
Eq. (|19p . This corresponds to a continuous single-particle 
density of states that grows quadratically. In Fig. [5l we 
show the pattern of complex zeros of the corresponding 
Z N (fi) for N = 50 and N = 100. Comparison with Fig.g] 
shows considerable difference: the number of zeros being 
much smaller for the case of continuous density of states. 
In the latter, /3 C is about 10 percent higher. In both, the 
estimated condensation temperaure T c increases approx- 
imately as N 1 / 3 . 



Fisher pointed out that there are complex zeros of the 
canonical partition function Zpj(0) on the complex /3 
plane. For a fixed N, the number of zeros on the com- 
plex plane, denoted by /3 r , is finite. Therefore Zm{P) 

may be expressed as a finite product Yl ~ ^ ms 

is unlike the Lee- Yang zeros whose validity was contin- 
gent on short-range interparticlc repulsion. At a phase 
transition, the complex Fisher zeros close in on the real (3 
axis. The Hclmholtz free energy for an A-particle system 

is given by F — — — lnZjv(/3). Since Zn(/3) is a sum of 
P 

exponential positive terms in (3, and is larger than unity 



IV. CONCLUDING REMARKS 

We have shown that the advent of a phase transition 
in a system is reflected in the pattern of the complex ze- 
ros of the partition function. Strictly speaking, a phase 
transition takes place only in the thermodynamic limit. 
But even for a finite system with relatively small number 
of particles, the pattern of complex zeros begin to close 
in on the real fugacity or inverse temperature axis. For 
the grand partition function Lee- Yang zeros, it is imper- 
ative to have short-range repulsion in the interparticlc 
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FIG. 5. The Fisher zeros for systems of 50 or 100 trapped 
bosons with a continuous energy distribution. 



interaction, whereas for the Fisher zeros of the canoni- 
cal partition function, this is not necessary. In the case 
of BEC, a signal of a phase transition is a peak in the 
heat capacity per particle on the real temperature axis, 
as shown in Fig. 3. This peak shows up nicely even for 
(N) = 50. In Fig. 4, the complex Fisher zeros for N — 50 
appear to close in on the real axis at the same tempera- 
ture. The Lee- Yang zeros in the van der Waal gas seem 
to close in towards the real z axis for T < T c . Finally, 
we note that even though there are no Lee- Yang zeros for 
the ideal Bose gas, the grand canonical partition function 
has a simple pole at z = 1, which was already noticed by 
Kastura.— 
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Appendix A: Fisher zeros 

The Fisher zeros are the zeros of the ./V-particle canon- 
ical partition function Zn(/3) in the complex (3 plane. 

Given Z (/3) = 1 and Z x (/3) = J2n=a e~^ En , we can 



obtain Z/v(/3) by the recursio: 

N 



,19 



1 N 



(Al) 



fc=i 



If Z\{j5) is the canonical partition function of a single 
particle in a three-dimensional harmonic oscillator well, 
then 



\n=0 



1 - e 



-Phuj/2 



-(3hw 



(A2) 



where e n = (n + l/2)hw. The calculation of Z^[fi) can 
be simplified by introducing 20 so that 



z N {y) - 



..3N/2 



n v =i(i-y) : 



■Pn(v). 



(A3) 



The recursion relation ()Alj) is reformulated as Po{y) 
Pi{y) = 1 and 



N 



E 

fe=i 



(i-y*) 3 



P N - k (y). (A4) 



The Pn(v) is a polynomial in y and when it is zero so is 
Zjq(y). In the case that e n = nhio rather than (n+ 1/2) Hlo 
we have 



Zi(y) = 



i 



(i-y) 2 



Z N (y) = 



1 



nf=i(i-y) 3 



P*(y), 



(A5) 

where the Pjv satisfy the same recursion (|A4|) . Thus the 
Fisher zeros will be the same irrespective of £rj = or 
huj/2. 

Since P/v(j/) is a polynomial the number of zeros is 
equal to its degree which increases rapidly with particle 
number. For example, for N = 50 there are 3495 zeros. 
We determine a subset to give a clear indication that the 
pattern of zeros pinches the positive real (3 axis. 

In order to numerically obtain the zeros we use the La- 
guerre method since this method "is guaranteed to con- 
verge to a (zero) from any starting point."— Once a zero 
is found we deflate the polynomial to obtain a second 
distinct zero and repeat. For larger N the zeros may be 
close together so we check the accuracy of the zero yo by 
considering a small circle centred on yo in the complex 
y plane and by ensuring that the curves Re(Pjv(2/)) = 
and Im(P/v(y)) = intersect inside the small disc. By 
using a radius of say 10~ 5 we have an estimate of the 
precision of the zero. Since we are interested in zeros 
that pinch the positive real (3 axis, we limit the variable 
y so that \y\ < 1 which results in Re(/3) > 0. 
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